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EFFICIENT  RECONSTRUCTION  OF  BLOCK-SPARSE  SIGNALS 


Joel  Goodman  * 

Naval  Research  Laboratory 
4555  Overlook  Ave.  SW 
Washington,  DC  20375 


ABSTRACT 

In  many  sparse  reconstruction  problems,  M  observations  are 
used  to  estimate  I<  components  in  an  N  dimensional  basis, 
where  N  >  M  ^  K.  The  exact  basis  vectors,  however, 
arc  not  known  a  priori  and  must  be  chosen  from  an  M  x  N 
matrix.  Such  underdetermined  problems  can  be  .solved  using 
an  £2  optimization  with  an  £1  penalty  on  the  sparsity  of  the 
solution.  There  are  practical  applications  in  which  multiple 
measurements  can  be  grouped  together,  so  that  K  x  P  data 
must  be  estimated  from  M  x  P  observations,  where  the  £1 
sparsity  penalty  is  taken  with  respect  to  the  vector  formed  us¬ 
ing  the  £2  norms  of  the  rows  of  the  data  matrix.  In  this  paper 
we  develop  a  computationally  efficient  block  partitioned  ho- 
motopy  method  for  reconstructing  K  x  P  data  from  M  x  P 
observations  using  a  grouped  sparsity  constraint,  and  compare 
its  perfonnance  to  other  block  reconstruction  algorithms. 

1.  INTRODUCTION 

Sparse  reconstruction  is  essential  to  compressed  sensing 
applications,  where  it  has  been  shown  that  it  is  possible  to  re¬ 
cover  a  /f-sparse  signal  of  length  N  from  0(I<  log  N)  com¬ 
pressive  measurements  1 1-3]. 

In  applications  where  measurements  are  obtained  from 
multiple  scn.sors  [4],  it  is  advantageous  to  group  the  sensor 
measurements  together  to  provide  robustness  against  impair¬ 
ments  such  as  noise  and  signal  fading.  The  formulation  for 
multichannel  measurements 

N 

luin  ||vec(r)  -  <I>vec(A')||^  +  /'  H  lle/’A'Ib  (1) 

i=l 

u.ses  an  £1  sparsity  constraint  vector  grouped  across  the  £2 
norms  of  the  rows  of  the  data  matrix,  where  <I>  G  jjAZ/'xwr 
is  the  dictionary  matrix,  the  matrix  Y  6  represents  M 

ob.servations  from  each  of  P  scn.sor  channels.  A'  6  is 

the  sparse  data  matrix,  vec(  )  is  an  operator  which  stacks  the 
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rows  of  a  matrix  to  fonn  a  column  vector  and  ej  6  is  a 
column  vector  of  Os  with  a  1  in  the  tth  position. 

To  solve  (1),  the  subgradients  are  formulated  with  respect 
to  blocks  (rows)  of  the  matrix  A'’,  atid  the  regularization  pa¬ 
rameter  fi  is  systematically  reduced,  tracking  the  solution 
A(/i),  until  a  new  £2  normed  row  of  the  solution  is  about  to 
turn  nonzero.  Reduction  in  /i  continues,  activating  more  and 
more  rows  of  A'(/r).  We  will  demonstrate  that  an  approx¬ 
imate  form  of  block  homotopy  outixirforms  greedy  block 
reconstmetion  (with  arbitrary  <h  in  (I))  for  roughly  the  same 
computational  complexity,  while  solving  (1)  for  all  sparsity 
levels  of  X.  The  rest  of  this  paper  is  organized  as  follows.  In 
Section  2,  we  extend  the  homotopy  technique  in  [5]  to  block 
partitioning.  In  Section  3,  we  demonstrate  the  performance 
of  block  homotopy  continuation  with  respect  to  other  block 
reconstmetion  algorithms,  and  in  Section  4  we  provide  a  brief 
summary. 

2.  BLOCK  PARTITIONED  HOMOTOPY 

We  return  to  the  original  optimization  problem  in  (1)  with 
a  slight  variation  in  notation, 

Hx)  ||y  -  <l>.rHi  +  /'  l|a--(^-)ll2,  (2) 

k=l 

where  x{k)  is  the  Ath  row  of  A',  and  y  e  and  x  6 

jjAf’xi  ,||.g  vec(Y)  and  vcc(A'),  respectively.  Complex  vec¬ 
tors  can  be  handled  in  the  form  (2)  by  grouping  the  real  and 
imaginary  components  together  after  observing  that  complex 
£2  norms  can  be  written  as  real  £2  norms  in  these  components. 

The  homotopy  approach  solves  the  optimization  problem 
based  on  (2)  for  large  /( and  then  tracks  the  solution  as  //  is  de¬ 
creased.  Initially,  /i  is  cho.scn  to  be  2||<I>^t/||2,  which  yields  an 
all-zero  solution.  Solutions  are  tracked  by  .solving  the  subgra¬ 
dient  c(|uation  fora  minimum.  Recall  that  the  .subdifferential 
at  X  of  a  convex  function  f{x)  is  the  .set  of  subgradient  vectors 
(  that  satisfy 

f(x)  -  /(x)  >  (x  -  (3) 

When  /(x)  is  differentiable  at  x,  the  subgradient  consists  of 
the  derivative  vector.  When  x  minimizes  f(x),  0  is  a  subgra- 
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client  of  f{x)  at  x.  The  derivative  of  ||x||2  is  x/||x||2  away 
from  zero.  With  /(•)  =  1|  •  II2  and  x  =  0,  we  observe  that 
||a-i|2  >  where  ||^||2  <  1. 

If  /(x)  is  the  sum  of  convex  functions  /  =  A 
each  of  which  has  subgradients  then  ^ 
gradient  of  /.  Using  this  observation,  the  subdifferential  of 
L{x)  contains  the  set 

{2($'^a>x  -  +  fi  [n"f  •  •  •  (4) 


where 


Uk 


dcf 


lkWI|2 

L  ■  ii^c-iu  <  1 


xik)  ^  0 
x{k)  =  0 


(5) 


Any  ^k  satisfying  the  constraints  can  be  used. 

Let  Cofr  represent  a  subset  of  the  positive  integers  less  than 
or  equal  to  N  such  that  k  e  Cou  implies  x{k)  =  0.  Let  Con 
represent  the  complement  of  Cqu  in  the  same  set  of  positive 
integers.  Without  loss  of  generality,  assume  that  the  compo¬ 
nents  are  ordered  so  that 


Thus,  letting  pk  =  Il-'i'lli  >  ll-'t-ib-  Similarly,  letting 

Pk  =  ||x(A-)||2/||ai|2,  we  have 

X]IWA-)1|2>IN|2.  (12) 

k 

Consequently,  the  nonnegative,  convex  function 

Liow  ='  111/  -  'f-'clli  +  f'lla-ib  (13) 

is  a  lower  bound  on  L(x).  If  we  show  Ljow  is  minimized  at 
X  =  0,  it  follows  that  x  =  0  also  minimizes  L(x).  Further¬ 
more,  this  minimum  is  unique  since  L|ow  is  strictly  convex 
except  potentially  on  rays  from  the  origin  where  l|.x||2  is  not 
strictly  convex.  However,  restricted  to  any  ray,  Liow  has  a 
unique  minimum  and,  hence,  x  =  0  is  the  unique  minimum 
ofL. 

To  minimize  Liow,  take  the  gradient  away  from  x  =  0  and 
take  the  inner  product  with  x,  getting 

(x’^VLiow)(a;)  =  2x’^(4>^4>x  -  ^’^y)  +  px'^ii,  (14) 


where  the  subscripts  refer  to  the  sets  Con  and  CoK-  Partition 
the  matrices  accordingly,  letting 


Then  zero  is  a  subgradient  if 


The  block  components  14.  with  k  e  Co,,  are  required  to  have 
unit  (2  norm,  even  if  the  corresponding  x(A)  vanishes.  This 
condition  holds  for  the  initial  nonzero  solution  and  will  hold 
throughout  the  homotopy.  We  can  unravel  the  matrix  equation 
to  get 


where  ii  =  x/||xl|2.  This  is  strictly  positive  if 

px'^u  =  /t||xH2  >  2|x'^<I'^(/|.  (15) 

Since  |x^<l>^(/|  <  l|xl|2||4>^i/||2,  we  have  a  strictly  positive 
inner  product  when  p  >  2||<I>^  J/II2.  In  this  case,  a  nonzero  x 
cannot  minimize  L\ow,  so  x  =  0  provides  the  minimum. 
Starting  with  x  =  0  and  p  =  2||<I>'^'i/||2  =  Hzofilb  = 
II2II2,  it  is  clear  that  (10)  is  solved  with  a  valid  t/on  = 

Ihu  lb  <  1-  la  (act,  this  is  the  case  as  long  as 

/!  >  max  ||2(A)||2-  (16) 

When  p  is  reduced  until  equality  holds  in  (16),  we  can  incor¬ 
porate  the  index  k  achieving  equality  into  the  .set  Con-  Etjua- 
tion  (10)  still  holds  since  Uk  =  z{k)/p.  Then  p  can  be  re¬ 
duced  further  using  (10).  The  sets  Co„  and  Coti  arc  modified 
in  the  following  manner. 

1.  If  reducing  p  causes  a  block  comimnent  14.  with  k  € 
Coir  to  achieve  unity  norm,  then  enter  the  k"'  block  into 
Con  and  reduce  Cotr  accordingly. 


•'^’011  —  ■'4  (2^0,,  /o/o,i)  (10) 

UqH  —  “(^olT  —  H  /i  2o„)  -f-  H  /I  tton- 


2.  If  reducing  p  cau.ses  the  k"‘  block  component  of  .2:0,1  to 
vanish,  then  add  the  k"‘  block  to  Cofr  anti  modify  Co„ 
accordingly. 


For  sufficiently  large  p,  (2)  has  the  all-zero  solution.  In 
fact,  for  real  pk,  observe  that 

min  IC'-  I  =  CTf  ,  ^ 

l|p|b=i'^  VEtlPil=i  / 

since  11-112  isaconvex  function  and  thus  assumes  its  maximum 
at  the  vertices  of  the  simplex  {p  :  pjt-  >  0  and  YlPk  =  (}- 


2.1.  A  Simplified  Approximaliun 

The  steps  evaluated  below  (for  the  case  when  the  block 
size  is  greater  than  unity)  assumes  Uonip)  is  piecewise  con¬ 
stant  between  branches.  Only  rule  I  applies  during  the  ho¬ 
motopy.  Wc  justify  this  by  noting  that  rale  2  involves  finding 
a  zeroed  block  component  of  Xo„  using  the  first  equations  of 
(10),  assuming  a  constant  u„„. 


The  new  nonzero  coinponenl  al  eacli  branch  is  deter¬ 
mined  by  (10);  once  activated,  this  block  component  remains 
fixed  tliroughout  the  algorithm.  Defining  to  be  a  matrix 
whose  orthononnal  columns  span  the  entries  associated  with 
the  block,  and 

r>T  1  —  1  r>T  4— 1 

O'  — ■  Zoff  B  2^011  j  ft  —  BA.  U()iu 

Ihe  norm  crossing  of  rule  1  is  expressed,  in  the  block  com¬ 
ponent  of  by  the  equality 

Wtlfh+O/kWl  =  II^^J(/'/^  +  rt)||2  = 

where 

ipT  ri  tiff  tnT/j 

au  =  El  a  fh-  =  El  fh 

This  quadratic  in  /( is  solved  for  each  block  component  using 

±  -  4|ky,||i(||/4|li  -  1) 

Tire  component  crossing  unity  norm  first  (largest  /<)  is  acti¬ 
vated.  Pseudocode  for  the  approximation  is  provided  in  Al¬ 
gorithm  1.  This  method  differs  from  another  approximation 
to  group  LASSO  [6]  in  its  solution  path  computation  as  well 
as  its  final  result. 


Algoriflllll  1  APPROXfMATEBLOCKHOMOTOPY 

Input:  <h  €  ij  €  w  e  Z 

k  «-  argmax^.  ||£*..$'^(/||2 

Ti  *”  W^x^'^yh 

Con  <-  {k}.  Con  <-  { j  e  Z|0  <  i  <  N,  i  5-^  k] 
while  IIj/  -  <I>.r||^  //||a;||i  >  £<  and  i  <  ?,na,  do 

A  <-  2<t>J„«I>on.  '24'Jn<I>ofr 

^OII  *-  ^o(T  <-  2<i>j,fi/ 

l/„„  <-  BLKNORMAf.fZn(i(2ofr-/I^/l"‘2on)-£:„„,l/) 

ft  <  .Z^on  B  A  Zou 
for  all  k  6  Cofr  do 

Hk  <-  argniaxn<,„<p  ||qa.  -f  =  > 

end  for 

k  <-  argniaxA  yk,  y  «-  /'if 

^on  *  ^011  U  (k),  C„n  <—  C„(T  \  {k}y  i  <—  )  +  1 

end  while 

rctum  Con 


3.  PERFORMANCE 

Monte  Carlo  simulations  were  iiin  comparing  the  de¬ 
tection  (dictionary  column  identification)  performance  of 
exact  and  approximate  block  pailitioncd  homolopy  process¬ 
ing,  as  well  as  S-OMP  [7].  The  signals  used  were  sparse 


in  a  Fourier  basis,  where  $  =  0lF,nu  6  cioox.soo^ 

€  C'’’™’'®™  representing  an  inverse  DFT  matrix,  and 
0  €  a  matrix  whose  elements  were  drawn  from  a 

zero  mean  unit  variance  Gaussian  distribution.  Monte  Carlo 
simulations  consisted  of  1000  runs  at  sparsity  levels  rang¬ 
ing  from  I  to  20  percent  in  1  percent  steps.  For  both  block 
partitioned  homotopy  and  S-OMP,  signals  were  blocked  in 
groups  of  P  =  5  (e.g.,  a  5-sensor  array).  Block  partitioned 
homotopy  processing  outperformed  S-OMP,  as  illustrated  in 
Fig.  1. 


Fig.  1.  Comparison  of  detection  performance  between  ap¬ 
proximate  block  partitioned  homotopy  continuation,  S-OMP, 
group  LASSO,  and  SOCP  where  the  signals  are  down  sam¬ 
pled  by  a  factor  of  5.  In  all  cases  the  SNR  was  30dB. 

We  also  compared  the  performance  of  approximate  block 
partitioned  homolopy  to  that  of  SOCP  using 

p 

min  V  ||.t(A:)||2  s.l.  ||y  -  <I>.'c||2  <  c,  (18) 

X  *• 


and  group  LASSO  with  the  block  shrinkage  function  (.soft  it¬ 
erative  thresholding) 


5,(:f(A-)) 


0 


if  Mk)h  >  ^ 
if  ||.T(A-)|l2  <  I  ' 

(19) 


replacing  the  shrinkage  function  given  in  [8].  In  the  cases 
simulated,  the  identification  iierformance  of  bltKk  partitioned 
homotopy  processing  was  virtually  identical  to  that  of  SOCP 
and  group  LASSO  for  the  discrete  sparsity  values  of  2.5%, 
5%,  10%,  15%  and  20%;  however,  there  was  a  iterformance 
advantage  to  block  partitioned  homotopy  processing  dis¬ 
cussed  in  the  computational  complexity  subsection.  Note 
that  group  CoSaMP  (9)  was  not  used  in  the  compari.son  as 
it  did  not  delect  signals  beyond  the  7%  sparsity  level.  This 
is  because  during  the  support  merger  phase  of  CoSaMP  3x 
as  many  columns  are  used  in  the  estimation  process  prior  to 
pruning. 


3.1.  Array  Processing  Application  Example 

We  simulated  a  CS  receiver  witli  P  =  4  aiUeniias  oper¬ 
ating  in  an  environment  with  RF  emitters,  including  hoppers, 
in  the  2-6  GHz  frequency  range.  The  4-antcnna  linear  ar¬ 
ray  was  spaced  in  2.5  cm  increments,  with  the  response  of 
each  antenna  constant  (unity  gain)  in  azimuth.  Compressive 
measurements  in  a  sensor-frequency  basis  were  obtained  via 
random  sampling,  with  4>  G  where  M  =  64  and 

N  =  512  (see  [4]  for  a  full  description  of  the  simulated  sce¬ 
nario).  As  shown  in  Fig.  2,  approximate  block  homotopy 


0{MNF^)  operations,  while  the  computation  of  re¬ 

quires  0{MNKP^)  operations  to  obtain  the  matrices  A  and 
B.  In  some  cases,  however,  there  may  be  sufficient  memory 
to  store  2<l>^4'  for  all  realiziitions  of  <I>. 

S-OMP’s  complexity  is  dominated  by  front  end  correla¬ 
tions  requiring  0{P^KMN)  operations  after  K  iterations, 
given  efficient  rank- 1  iqxlates  are  used  to  compute  the  projec¬ 
tor  used  in  reconstruction.  SOCP  has  complexity  0[{KN)^) 
(1 1 J,  and  in  general,  group  LASSO  approaches  to  solving  (1) 
require  iteratively  searching  for  a  /r  or  e  to  achieve  the  desired 
level  of  sparsity. 


SNR(dB) 

Fig.  2.  Comparison  of  detection  performance  between  ap¬ 
proximate  block  partitioned  homotopy  continuation  and  S- 
OMP  as  a  function  of  SNR  in  an  array  processing  application. 
Note  that  the  signals  on  each  of  the  4  channels  is  randomly 
down-sampled  by  a  factor  of  8. 

detection  performance  was  noticeably  better  than  that  of  S- 
OMP  at  all  SNRs  and  sparsity  levels  measured. 

3.2.  Computational  Complexity 


4.  SUMMARY 

In  this  paper  we  present  block  partitioned  homotopy  pro¬ 
cessing  for  multichannel  sparse  signal  reconstruction.  We  ex¬ 
tended  the  homotopy  continuation  [5]  to  basis  pursuit  opti¬ 
mization  with  a  block  partitioned  sparsity  constraint  that  con¬ 
sists  of  the  C2  norm  of  each  row  of  the  signal  matrix  combined 
in  an  fashion.  Block  partitioned  homotopy  processing  out- 
perfonned  S-OMP  in  identifying  the  columns  vectors  of  the 
dictionary  spanned  by  signal,  and  did  so  with  roughly  the 
same  computational  complexity  as  S-OMP.  Block  partitioned 
homotopy  processing  performance  w.is  on  par  with  SOCP,  but 
with  a  significantly  lower  computational  complexity. 
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